C 34567============================================================70
C     Divide the compute domain into NNPART parts
C     Solve the control equation system in each part
      SUBROUTINE DCMD_MAIN(OPT, NNODE, L, W, T, P, V, 
     &                     TW, JH, JM, TPC, EFF, IFLAG)
C     Invoke global declaration
      USE DCMD_CMOD
      
      IMPLICIT NONE
      
C     OPT - (input) options of concurrent or countercurrent case
C     NNODE - (input) node number to restore local properties
C     L - (input) length of membrane
C     W - (input/output) 2D array to restore local mass flowrates
C     T - (input/output) 2D array to restore local temperatures
C     P - (input/output) 2D array to restore local prssures
C     V - (input/output) 2D array to restore local velocity
C     TW - (output) 2D array to restore local wall temperatures
C     JH - (output) 1D array to restore local heat transfer flux
C     JM - (output) 1D array to restore local mass transfer flux
C     TPC - (output) 1D arary to restore local temeprature 
C                       polarization coefficient
C     EFF - (output) 1D array to restore local thermal efficiency
C     IFLAG - (output) integer flag of running status
C             IFLAG = -1 means that fail to get the initial value
C             IFLAG = 0 means that the iteration was not run
C             IFLAG = (any other integer) means that the number of
C                     iteration
      INTEGER*4 OPT, NNODE, IFLAG
      REAL*8 L
      REAL*8 W(2,NNODE), T(2,NNODE), P(2,NNODE), V(2,NNODE)
      REAL*8 TW(2,NNODE)
      REAL*8 JH(NNODE), JM(NNODE), TPC(NNODE), EFF(NNODE)
C     Intermediate variables
      INTEGER*4 I, ICOUNT, ISIDE, NMAX, NPART
      REAL*8 TRELERR, TRELERR0, RBT(2), LBT(2), TSTEP, DRBT
      REAL*8 WRELERR, WRELERR0, RBW(2), LBW(2), WSTEP, DRBW
      REAL*8 WINLOC(2), TINLOC(2), PINLOC(2), VINLOC(2),
     &       WOUTLOC(2), TOUTLOC(2), POUTLOC(2), VOUTLOC(2)
      REAL*8 LENLOC
      REAL*8 DAMP
C     Function declaration
      REAL*8 GradDP
        
      NMAX = 1000
      NPART = NNODE-1
      
      LBT(1) = T(1,1)
      LBT(2) = T(2,1)
      RBT(1) = T(1,NNODE)
      RBT(2) = T(2,NNODE)
      TSTEP = 0.D0
      
      LBW(1) = W(1,1)
      LBW(2) = W(2,1)
      RBW(1) = W(1,NNODE)
      RBW(2) = W(1,NNODE)
      WSTEP = 0.D0
      
      DAMP = 0.5D0
C     A extreme large initial value of RELERR0 is set to ensure that
C     the convergent judge is valid in the first iternation.
      TRELERR0 = 1.D5
      WRELERR0 = 1.D5
      
      LENLOC = L/NPART
      
      SELECT CASE(OPT)
C     Concurrent case
      CASE(1)
        DO I = 1, NNODE-1
          DO ISIDE = 1, 2
            WINLOC(ISIDE) = W(ISIDE,I)
            TINLOC(ISIDE) = T(ISIDE,I)
            PINLOC(ISIDE) = P(ISIDE,I)
            VINLOC(ISIDE) = V(ISIDE,I)
            COM_WLOC(ISIDE) = WINLOC(ISIDE)
            COM_TLOC(ISIDE) = TINLOC(ISIDE)
            COM_PLOC(ISIDE) = PINLOC(ISIDE)
            COM_VLOC(ISIDE) = VINLOC(ISIDE)
          END DO
          CALL DCMD_SOLV(OPT, LENLOC, WINLOC, TINLOC, PINLOC, VINLOC,
     &                   WOUTLOC, TOUTLOC, POUTLOC, VOUTLOC, 
     &                   TW(:,I), JH(I), JM(I), TPC(I), EFF(I))
          DO ISIDE = 1, 2
            W(ISIDE,I+1) = WOUTLOC(ISIDE)
            T(ISIDE,I+1) = TOUTLOC(ISIDE)
            P(ISIDE,I+1) = POUTLOC(ISIDE)
            V(ISIDE,I+1) = VOUTLOC(ISIDE)
          END DO
        END DO

C     Countercurrent case
C     Left boundary of compute domain is the inlet of tube side
C     Right boundary of compute domain is the inlet of shell side 
      CASE(2)
        DO ICOUNT = 1, NMAX
          IFLAG = 0
          DO I = 1, NNODE
            DO ISIDE = 1, 2
              WINLOC(ISIDE) = W(ISIDE,I)
              TINLOC(ISIDE) = T(ISIDE,I)
              PINLOC(ISIDE) = P(ISIDE,I)
              VINLOC(ISIDE) = V(ISIDE,I)
              COM_WLOC(ISIDE) = WINLOC(ISIDE)
              COM_TLOC(ISIDE) = TINLOC(ISIDE)
              COM_PLOC(ISIDE) = PINLOC(ISIDE)
              COM_VLOC(ISIDE) = VINLOC(ISIDE)
            END DO
            CALL DCMD_SOLV(OPT, LENLOC, WINLOC, TINLOC, PINLOC, VINLOC,
     &                     WOUTLOC, TOUTLOC, POUTLOC, VOUTLOC, 
     &                     TW(:,I), JH(I), JM(I), TPC(I), EFF(I))
C     New tube-side values are produced in the right side of given
C     values in each iteration.
            ISIDE = 1
            IF (I .LT. NNODE) THEN
              W(ISIDE,I+1) = WOUTLOC(ISIDE)
              T(ISIDE,I+1) = TOUTLOC(ISIDE)
              P(ISIDE,I+1) = POUTLOC(ISIDE)
              V(ISIDE,I+1) = VOUTLOC(ISIDE)
            END IF
C     New shell-side values are produced in the left side of given
C     values in each iteration.
            ISIDE = 2
            IF (I .GT. 1) THEN
              W(ISIDE,I-1) = WOUTLOC(ISIDE)
              T(ISIDE,I-1) = TOUTLOC(ISIDE)
              P(ISIDE,I-1) = POUTLOC(ISIDE)
              V(ISIDE,I-1) = VOUTLOC(ISIDE)
            END IF
          END DO
C     Outlet pressure calculation
          P(1,NNODE) = COM_PIN(1)
     &    +COM_LEN*GradDP(COM_DEQ(1), COM_RHO(1), COM_VIN(1), COM_MU(1))
          P(2,1) = COM_PIN(2)+COM_LEN*
     &             GradDP(COM_DEQ(2), COM_RHO(2), COM_VIN(2), COM_MU(2))
          
C     The right boundary properties shall be varied in an enough small
C     range, when it is convergent.
          TSTEP = T(1,NNODE)-RBT(1)
          WSTEP = W(1,NNODE)-RBW(1)
          TRELERR = DABS(TSTEP)/RBT(1)
          WRELERR = DABS(WSTEP)/RBW(1)
          IF (TRELERR .LT. 1.D-5) THEN
            IFLAG = IFLAG+1
          ELSE
            RBT(1) = T(1,NNODE)
            IFLAG = 0
          END IF
          IF (WRELERR .LT. 1.D-5) THEN
            IFLAG = IFLAG+1
          ELSE
            RBW(1) = W(1,NNODE)
            IFLAG = 0
          END IF
C     The left boundary properties shall be varied in an enough small
C     range, when it is convergent.
          TSTEP = T(2,1)-LBT(2)
          WSTEP = W(2,1)-LBW(2)
          TRELERR = DABS(TSTEP)/LBT(2)
          WRELERR = DABS(WSTEP)/LBW(2)
          IF (TRELERR .LT. 1.D-5) THEN
            IFLAG = IFLAG+1
          ELSE
            LBT(2) = T(2,1)
            IFLAG = 0
          END IF
          IF (WRELERR .LT. 1.D-5) THEN
            IFLAG = IFLAG+1
          ELSE
            LBW(2) = W(2,1)
            IFLAG = 0
          END IF
          
          IF (IFLAG .EQ. 4) THEN
            IFLAG = ICOUNT
            EXIT
          END IF
          
          IF (ICOUNT .EQ. NMAX) THEN
            WRITE(*, *) 'Error: Iteration cannot be convergent!'
          END IF
          

!          TSTEP = T(2,1)-LBT(2)
!          TRELERR = DABS(TSTEP)/LBT(2)
!C         Judge the convergence in solution
!C         If the relative error is larger and larger,
!C         it indicates divergence
!          IF (TRELERR .GE. TRELERR0) THEN
!            IFLAG = -1
!            STOP
!          ELSE
!            TRELERR0 = TRELERR
!            LBT(2) = T(2,1)
!          END IF
!          
!C         The new trying value of T(2,1) was determined by the value 
!C         of TSTEP. 
!          IF (TRELERR .GE. 1.D-5) THEN
!C           DRBT - the range of T(2,1) changing
!            DRBT = T(1,1)-T(2,1)
!C           Ensure that the new trying value of T(2,1) is not out of
!C           the range of DRBT
!            IF (DABS(DRBT) .LT. DABS(TSTEP)) TSTEP = DSIGN(DRBT, TSTEP)
!            T(2,1) = T(2,1)-TSTEP*DAMP
!            IFLAG = ICOUNT
!          ELSE
!            IFLAG = ICOUNT
!            EXIT
!          END IF
!C
!          WSTEP = W(2,NNODE)-LBW(2)
!          WRELERR = DABS(WSTEP)/LBW(2)
!C         Judge the convergence in solution
!C         If the relative error is larger and larger,
!C         it indicates divergence
!          IF (WRELERR .GE. WRELERR0) THEN
!            IFLAG = -1
!            STOP
!          ELSE
!            WRELERR0 = WRELERR
!          END IF
!C         The new trying value of W(2,1) was determined by the value 
!C         of WSTEP. 
!          IF (WRELERR .GE. 1.D-5) THEN
!C           DRBW - the range of W(2,1) changing
!            DRBW = W(2,1)-W(2,NNODE)
!C           Ensure that the new trying value of T(2,1) is not out of
!C           the range of DRBT
!            IF (DABS(DRBW) .LT. DABS(WSTEP)) WSTEP = DSIGN(DRBW, WSTEP)
!            W(2,1) = W(2,1)-WSTEP*DAMP
!            IFLAG = ICOUNT
!          ELSE
!            IFLAG = ICOUNT
!            EXIT
!          END IF

        END DO

      END SELECT
      
      RETURN
      
      END SUBROUTINE
